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Abstract 


BLISS (Bi-Level Integrated System Synthesis) is a decomposition optimization 
method for engineering systems. The method is characterized by the separate 
optimization of a relatively few system-level variables and the optimization of potentially 
numerous local variables. Subsystem optimizations are autonomous and may be 
conducted concurrently (i.e. on a multiple processor computer). In previous versions of 
BLISS, optimum sensitivity analysis and system sensitivity data were used to link the 
subsystem optimization data to the system optimization. The current work replaces both 
the optimum sensitivity analysis and the system sensitivity equations by the quadratic 
response surface representations using subsystem optimization results. 

The response surface methodology for BLISS achieves the desired improvements 
while retaining key attributes of previous versions of BLISS: the autonomy of the black 
box optimizations and the clear separation of the system variables from the potentially 
numerous local variables. The response surface formulation of BLISS was successfully 
demonstrated on a simplified conceptual design of a supersonic business jet. 

In addition to changes in the overall optimization methods of BLISS, subsystem 
fidelity was enhanced, accompanied by the necessary modifications to the data flow 
between subsystem analyses. Documentation of these modifications that have not as yet 
been tested is included as a reference for future research. 
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1. Introduction 


The design of a complex engineering system (e.g. an aircraft) involving a large 
number of subsystems, modules, or black boxes (BB’s) inherently involves a large 
number of design variables and constraints. The optimization problem can quickly 
become too large to manage efficiently, and the solution process can become extremely 
expensive in a computational sense [1], Decomposition of the problem into more 
manageable subtasks enables an efficient distribution of work across disciplines. In such 
decomposition, the design variables and constraints local to a particular module are 
separated from those that affect the system as a whole. The separation and distribution of 
work promotes a diverse grouping of human and computer resources, thus conforming to 
current trends in parallel processing technologies and concurrent engineering. 

The original formulation of Bi-Level Integrated System Synthesis, BLISS 
(re (erred to as BLISS-98), introduced in [2] and documented in detail in [3], involved the 
conceptual design of a business jet, optimized for maximum range. During this design 
analysis, simple empirical and analytical calculations were performed for each of the BB 
analyses. In [4], the test case was the same supersonic business jet but the fidelity of the 
analyses was improved. Both previous BLISS formulations share the distinguishing 
feature that the system objective was range while that of local optimization was a 
composite function made up of the sum of BB outputs, each weighted by the system 
sensitivity derivatives for that term. Therefore, the optimum sensitivity derivatives are 
used as the coupling mechanism between system and local optimizations. 

The current version of BLISS (called BLISS-RS herein to denote response surface 
methodology) is conceptually similar to previous versions. However, in BLISS-RS the 
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coupling between system and local optimization occurs via quadratic response surfaces 
instead of an optimum sensitivity analysis. The new framework of BLISS was 
successfully tested using BB’s inherited from previous versions of BLISS. 

2. BLISS-98 System Architecture 

Regardless of the particular decomposition approach, the introduction of the 
BLISS algorithm begins with formulation of the problem without decomposition. A 
modular system like BLISS typically optimizes three different types of design variables. 
First, the system level variables Z affect at least two of the BB’s. Secondly, the local 
variables X are specific to a particular BB. Finally, the coupling variables Y A are output 
from a BB, while Y* are coupling variables input to a BB. A statement combining 
analysis and optimization of a generic system can be written as: 

Find: V Equation I 

Where V represents the collection of all design variables in 
the space {Z | X | Y* | Y A } 

Minimize: F(V) 

Satisfy: (g(V)} < 0, for each BB 

{h(V)}=0 
{c(V)} = Y* - Y A = 0 
{VL < X < VU} . 

In equation 1 , the inequalities {g} represent the behavior constraints local to a BB, and 
the equalities {h} correspond to the solution of the BB governing equations (BB inner 
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analysis). The coupling equalitites {c} describe the condition that the output coming 
from one BB has to match the input to another BB, for a particular coupling variable Y. 

Solution of the all-in-one formulation shown above can be very arduous due to 
the possibly large number of local variables. For this reason, it is beneficial to 
decompose the problem based on the types of design variables. BLISS is such a 
decomposition method in that the problem is divided into two steps: local optimization 
within each BB and the system optimization. 

The original formulation of BLISS relied predominately on gradients to guide the 
search toward an optimum. It depended on system and sensitivity analyses, local 
optimizations inside the BB’s, and the system optimization. Each cycle through the 
BLISS procedure improved the design in two steps; first by optimizing each BB for the 
local design variables X while holding Z constant; and next, a system-level optimization 
that treats Z as the design variables [3j,[4j. Figure 1 depicts the BLISS-98 process. 
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Figure 1 . BLISS-98 Cycle 


After initialization, the BLISS-98 performs a system analysis and sensitivity 
analysis in which Y and the derivatives of Y with respect to Z and X are computed. A 
linear approximation to the system objective (an element of Y namely Range output 
from the performance BB) as a function of Z and X is established using the above 
derivatives. That approximation is adopted as the objective function in BB optimizations 
that follow. In each BB, the Z and Y variables are frozen and an improvement in the 
objective function is achieved by the local optimizations that use local X separately in 
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each module. This is followed by computation of the derivatives of the optimum X with 
respect to the parameters Z and Y. 

The second step seeks improvement for the system-level variables Z and is linked 
to the first step by the derivatives of Xopt with respect to parameters Z and Y. The 
derivatives are used to extrapolate each subdomain optimum as a function of Z and Y. 
The functional relation Y=Y (Z) is approximated by extrapolation based on the system 
sensitivity analysis. These steps alternate until convergence. Note that the output of step 
1 is an optimum change in the local design variables, AXopt, in the presence of constant 
Z, and the output of step 2 is an optimum change in system design variables, AZopt. 

BLISS-98 was successfully implemented for the test case of a design of a 
supersonic business jet. A detailed description of the results obtained can be found in [4], 
However, its dependence on the system analysis and derivative information proved to be 
computationally costly. 

2.1. BLISS-98 System Variable Flow 

Figure 2 show's the general flow' of data between the various modules of BLISS- 
98. The system variables (those which affect more than one black box directly) are 
shown in the dashed boxes. Variables output from a BB are specified using ( A ), and 
those input to a BB are designated by (*). 
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System Variables Z 


Svstem Variables Z 


System Variables Z 



Figure 2. BLISS-98 System and Coupling Variable Flow 


2.2. BLISS-98 Subsystem Analyses 

The subanalysis modules included in BLISS-98, as well as the general data flow 
among these modules, were used as a benchmark for BLISS-RS. Therefore, a brief 
discussion of the analysis involved in each module is presented here. 

Structures BB 

Figure 3 is a schematic representation of the structures black box, showing system 
variable flow and intermediate operations. The structures BB takes in all Z-level 
variables excluding atmospheric parameters, along with Lift from the aerodynamics BB 
and engine weight from the power BB. 
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Figure 3. System Variable Flow for BLISS-98 Structures Black Box 


The aerodynamic loads are generated within the structures module in a pre- 
processor to generate the appropriate input into the structural analysis package. The lift 
loads on the wing are determined via the Shrenk approximation: a spanwise average of an 
elliptical lift distribution and a trapezoidal distribution that reflects the w'ing chord taper, 
as shown in Figure 4 [5]. This load distribution is based on half of the total aircraft 
weight. Lifting effects of the fuselage and tail are ignored. 
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Lift Distribution 



Figure 4. BLISS 1998 EL APS Lift Distribution [4] 


The structures BB uses as its primary analysis tool the Fortran-77 based 
Equivalent Laminated Plate Solution (ELAPS), described in [6-9]. ELAPS analyzes 
trapezoidal sections, or elements, that represent whole lifting surfaces (plate segments) or 
fuselages (shell segments). For each segment, a Ritz-based method is used to minimize 
strain energy, yielding polynomial equations for static deflections and internal stresses. 

Although the accuracy of ELAPS has been shown to be slightly below that of 
finite element codes, many beneficial aspects lend it nicely to multidisciplinary 
optimization. The main advantage lies in the fact that, compared to finite element 
codes, ELAPS generates far fewer degrees of freedom, thus reducing the required 
computational time and expense [10]. The input data is also much simpler and faster to 
develop. The model can consist of one segment representing an entire simple wing, or 
thousands of segments of a partitioned wing. Therefore, the adaptability of the model 
makes ELAPS an attractive option for all phases of the design process. 

Within the BLISS framework, stresses are analyzed along a three-segment 
wingbox. Each wingbox consists of the top and bottom sandwich panels of different 
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thicknesses and sandwich webs that are identical in the front and rear of the wingbox. 
The front spar of the wing box is located at 1 0% of the chord length and the rear spar lies 
at 70% of the chord length. Figure 5 depicts the configuration of the ELAPS model used 
by BLISS. 



The top and bottom panels as well as the webs have the thickness of the sandwich 
face sheets (tj) and the sandwich caliper thickness (%) as design variables, depicted in 
Figure 6. ELAPS models such a built-up structure by representing each face and the core 
as separate elements linked in a common coordinate grid. 



Figure 6. Wingbox Cross-section! [4] 
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The stresses calculated by ELAPS are used in formulation of a set of local 
constraints. The displacements are combined with the wingbox model to yield the value 
labeled “twist” in Figure 3. “Twist” is formulated as a change in the total lift of the wing 
divided by dynamic pressure, q, having units of ft 2 . Thus, it can be viewed as a change in 
effective wing area due to changes in the local chord angle of attack induced by the wing 
twist and bending. Finally, the total structural weight is calculated using the structures 
BB local variables. This is combined with the engine weight to estimate the total weight 
of the aircraft. Fuel weight is estimated by empirical relationships. 


Aerodynamics BB 



Figure 7 is a schematic representation of the aerodynamics black box, showing 
system variable flow and intermediate operations. The aerodynamics BB takes in all Z- 
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level variables, along with total weight and twist from the structures BB, and engine 
scaling factor from the power BB. 



Figure 7. System Variable Flow for BLISS-98 Aerodynamics Black Box 


The left side of the chart shows the inputs used by the AWAVE preprocessor. 
A WAVE is a simplified version of the Harris far-field wave drag program. The AWAVE 
Fortran code computes wave drag on the basis of the aircraft’s cross-sectional distribution 
along the centerline, hence it requires data about the entire configuration geometry to 
enable the area ruling [4]. 

A baseline drag polar is modified using ESF to account for wave drag of the 
engine nacelle. This is combined with the wave drag from AWAVE and drag due to lift 
to yield the wing’s drag coefficient Cowing- The drag coefficient for the horizontal tail, 
Coht, is calculated from an. empirical relation with Cowing- This leads to the value of lift to 
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drag ratio, L/D. Since BLISS models only the cruise flight regime, the total weight input 
set equal to total lift, L. Total drag is easily found dividing total lift by L/D. 

Power BB 

The power black box is a relatively simple in form. The drag input from the 
aerodynamics BB is set equal to thrust for the cruise condition. Engine weight is found 
through an empirical formula involving drag and ESF. SEC is determined using an 
engine deck approximation model (quadratic) that takes in values for thrust and throttle 
setting. 



Figure 8. System Variable Flow for BLISS-98 Power Black Box 
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Performance BB 


The simplest black box is the performance, or Range BB. It is made up entirely 
of the Breguet Range equation. 



Figure 9. System Variable Flow for BLISS-98 Range Black Box 


3. BLISS-RS 

The motivation for the current research is a direct consequence of the 
aforementioned shortcomings of BLISS-98. At the same time, the primary attributes of 
BLISS (clear separation of system and local design variables and the autonomy of BB’s) 
need to be preserved. Reformulating BLISS to incorporate response surfaces is a logical 
choice for addressing these problems. Since RS methods efficiently explore the entire 
design space, the optimization procedure is less likely to terminate at a local optimum. 
Also, the expense associated with calculating the local and system sensitivity derivatives 
is eliminated in the present method. However, this reduction in computational labor is 
partially offset by the labor to generate a number of response surfaces for each BB [11]. 
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Figure 10 shows the general BLISS-RS procedure. The first operation involves 
initialization of the system level variables Z, local variables X, and weighting factors w. 
Next, the coupling variables are initialized, either through a system analysis or an 
educated guess. After all variables are initialized, initial upper and lower bounds (LBo 
and UBq) are selected for all variables. A DOE pattern is then used to create a dispersion 
of inputs to each BB. The BB’s are optimized locally for each input, and then response 
surfaces are fitted through each optimal BB output. This collection of response surfaces 
is used for system optimization, allowing for rapid extraction of coupling data. If the 
system successfully converges, a final analysis is performed to retrieve local variables X. 
If convergence has not been achieved, the response surfaces are adjusted, and the system 
optimization process is repeated. 



Figure 10. BLISS-RS Flowchart 
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3.1. Local Optimization 


Formulation of a local BB objective function can be accomplished by observing 
the data flow between the various modules. In Figure 2, the system objective is the 
range, which is output from the performance BB. The range is directly related to those 
quantities entering the BB - namely total weight, fuel weight, specific fuel consumption, 
and lift-to-drag ratio. To illustrate an indirect influence of a BB output to the range, 
consider the wing twist. The twist information is sent to the aerodynamics module, 
where it is used to modify the wing shape and thus alter the lift and drag that the 
aerodynamics BB sends to the performance BB. 

Naturally, the structures BB should be optimized such that the range is increased 
as much as possible. However, some of the outputs of the structures BB influence 
performance directly, while others have an indirect effect. Therefore, it is not clear how 
much relative importance needs to be given to each output. 

This problem is addressed by assign ing weight factors w to each of the outputs of 
a BB, then optimizing for a composite objective function. These weight factors are 
analogous to the system sensitivity derivatives used in previous versions of BLISS. The 
composite objective function is a weighted sum of the contributions of each output. For 
example, the composite objective for the structures BB would be the sum of total weight, 
fuel weight, and wing twist, each multiplied by a weight factor. Since the actual values 
of w are unknown, a RS is formed in the space of (Z | Y* | w} and the task of 
determining the values of w is now left up the system optimizer. 

Given vectors of system variables, coupling variables, and weight factors, a black 
box can be optimized for a composite objective function, F with the satisfaction of all 
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local constraints. The result is a point F in the space {Z | Y* | w} . It represents the best 
contribution to the “synthetic” cost function the BB can make given Z, Y, and local 
constraints. The term “synthetic” is used because the system objective is not considered 
at this point; instead, the function that is minimized is made up of multiple terms that 
represent all outputs of the BB. Each of these terms includes a weighting factor w that 
dictates the relative importance of each response. The formation of the BB response 
surfaces (i.e. the BB optimizations) can be stated formally for BB as: 


Given: 

Z, Y*, w 

Find: 

X local to BB 

Minimize: 

ii 

M 

^ > 

Satisfy: 

{h} “ 0 


{g} * o 


(XL < X < XU} 


Where FT is the BB contribution to the system objective. The constraints {h} = 0 in 
equation 2 represent the analyses (FLAPS, for example). The {g} constraints could stand 
for physical limitations (material properties, dimensional tolerances, etc.). 

Using a DOE-based point placement method, Z, Y*, and w are varied to yield a 
locally optimized point in the space (Z | Y* | w} subject to local constraints. Many such 
points are generated to produce a cloud of points through which a quadratic response 
surface can be fit. A “surrogate” of each BB is obtained, and this approximation of each 
BB analysis can then be used in a system-level optimization. 
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3.2. System Optimization 


Having obtained the quadratic representations of each BB, the next BLISS action 
involves reading values from the quadratic RS’s to achieve an improvement in the system 
objective while satisfying the coupling constraints, (c}=0. The coupling constraints are 
easily recognized by the fact that the output BB; must be equal to the corresponding input 
to BBj. For example, the total lift generated in the aero BB must equal that passed along 
as an input to the structures BB. During system optimization, a coupling variable Y A is 
not directly sent from BBj as an input to BBj. Rather, the system optimizer proposes an 
input to BBj, Y*, and a constraint is introduced in the form Y A - Y* = 0. In Figure 11, 
the coupling constraints are shown graphically as circles where data flow channels meet. 



Figure 11. BLISS-RS System Optimizer 
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The system optimization can be formally stated as: 

Given: f>j for each BB Equation 3 

Find: Z,Y*, w 

Maximize: F = Range 

Satisfy: {c} = 0 

{g! =£ 0 

{XL < X < XU} 

In the above expression, (5j are the regression coefficients for the response surface 

Lh 

describing the j output of a particular BB. 

4. Response Surface Methodology 

The underlying goal of many types of experimentation is to correlate an output 
response to a set of factors of interest to the researcher. These relationships can be 
attained by constructing a model that describes the response over various values of the 
factors of interest. These response surfaces can be generated as a graphical means of 
displaying such relationships. Response surfaces can be analyzed to determine optimum 
combinations of input factors, or to explore relevant tradeoffs when multiple responses 
are involved. Figure 12 demonstrates such features. 
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Figure 12. Generic Example of a Two-inpul Response Surface 
With regard to design of a complex engineering system, the response surface 
methodology (RSM) allows the design space to be efficiently explored to determine the 
values of the design variables that optimize performance characteristics subject to system 
constraints [12]. RSM is used to obtain the mathematical models that approximate the 
functional relationships between performance characteristics and design variables [13]. 
Various design of experiments techniques, such as the central composite design (CCD) 
and the Box-Behnken design, are used to sample the design space efficiently [14]. With 
these approaches, experiments are performed at statistically selected locations in the 
design space. The resulting data is then used to construct response surface models 
through least squares regression. 
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However, these, like most other experimental designs were formulated with 
physical experiments in mind where the measurement variance of the response is the 
main concern. This stems from the fact that the output response of a physical system 
usually exhibits some degree of variability with the experiment repeated using the same 
inputs. Conversely, computer experiments (commonly used in engineering design) are 
deterministic, that is to say that there is no measurement error and no variability in 
response outputs given multiple runs of the same test using the same inputs [13], 
Therefore, experimental designs constructed to minimize variability of measurements 
may not be the best choice for computer experiments. 

4.1. Construction of Empirical Models 
4.1.1. Linear Regression 

In most applications of RSM, it is necessary to develop an approximation model 
to the true response surface. Approximation is necessary since, for most cases, the 
underlying function that drives the response is an unknown physical mechanism [12]. 
Multiple regression is used to generate an empirical model. 

A multiple linear regression model with n independent variables takes the form 

n 

y = 1 3 0 + ^ /l.y. + £ Equation 4 

i=i 

where y represents the response value, x,- are the independent variables ( predictor 
variables or regressors), pi are unknown partial regression coefficients, and e is the 
model error. This is a linear model because it is a linear function of the unknown 


20 


parameters fio, fii t P 2 , .... fin- This model describes a plane for two independent variables 
and a hyperplane in higher dimensions for three or more independent variables. 

Aany regression model that is a linear function of the regression coefficients is a 
linear regression model, regardless of the shape of the surface that it generates [12], For 
example, interaction terms could be added to equation 3 (n- 2) to give 

y = fi 0 + jSjjq + fi,x^ + fi 12 X 1 X 2 + £. Equation 5 

But if we let Xj=xiX 2 and fii~ fin, then equation 5 becomes a standard multiple linear 
regression model. Higher order models can be generated using similar techniques. 

The method of least squares is typically employed to evaluate the fi’s in equation 
5. If the response is observed m times, yj, yj, y m , where m>n, equation 5 can be 
rewritten as 


y i -Po + 2,P> x ij + s i> i=J - 2 m 

j - 1 

In matrix notation, this can be written as 

{ y}=[X]{fi} + {e } 

where 



Til 

I x n x 12 ■ 

' X ln 

ii 

y ‘ • 1^1 - 

1 x 21 x 22 

'• X 2n 


y m \ 

1 X ml X m2 ■ 

' X mn. 



A 


f i 

w- 

fi 2 

, and {e} = 

e 2 


A, 


£ m . 


Equation 6 


Equation 7 
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where { v} is an m-by- 1 vector of observations, [X] is an m-by-n matrix of the independent 
variables, {|3} is an n-by-1 vector of regression coefficients, and {e} is an m-by-1 vector 
of random errors. 

The least squares function is given by 

m 

{L} = Ixf ={e}’{e} = (M-[^]{P})’({f}-m{p}) 

i=l 

= {y}’{y} -mm y}-{yrmm+mxr[xm 

= M’{y}-2{p}’[X]’{y}+ {p nxnxm Equations 

Tire components of the random error vector are minimized by 

1^7 =-2[x]’{y}+2[mm} = 0 

which , after simplification, becomes 

[X]’[X]{b}=[X]’{y}. Equation 9 

Tire vector {6} represents the vector of least squares estimators of {P} that 
minimizes {L}. Finally, {6} is isolated by multiplying through by the inverse of [X]’[X], 
Tlris gives 

{b} = ([X]’[ X\) A [X]’{y] . Equation 10 

Tire fitted regression model is 

{y} = [X]{h}. Equation 11 

In indicial notation, Equation 1 1 is 

ft 

y i =b 0 +^ l b j X ij , i=l,2,...,m Equation 12 
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4.1.2. Second order models 


The techniques described above are easily adapted to approximate a second order 
model of the form 

n n 

y = Ao + ^ A 2 Pu x f + 22 P‘i X i X J +£ - Equation 13 

i-1 i-1 i< j 

In this case, the only difference lies in the formation of the [X] matrix. The 
minimum number of response observations needed to fit a second order model is given 
by 

(n 2 + 3n + 2) 

NS — Equation 14 

where n is the number of design variables. Put another way, the length of y must be at 
least (n 2 +3n+2)/2. The first column of [X] corresponds to the fio term in all observations, 
so all entries in this column are set to unity, file columns 2 through (n+ J) correspond to 
the linear terms. Columns (n+2) through (2n+ 1) correspond to the pure quadratic terms 
while the interaction terms are located from columns (2n+2) to ((n' : +3n+2)/2). Once the 
[X] matrix is formed, the vector {b} is found in the same manner as above, and a 
quadratic response surface is easily attained. 

The effect that dimensionality has on the total computational labor involved in 
generating response surfaces. Since the number of points needed to fit a quadratic 
increases with the number of design variables according to equation 14, efforts should be 
made to condense the number of design variables whenever possible. For a numerical 
example, consider a wing load distribution that is defined by forces at 10,000 discrete 
grid points on the wing surface in 2D coordinates (u,v). The approximate force p at a 
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point located at (u,v) is given by p = p(q,u,v). However, if p is formulated as a quadratic 
polynomial in q then the number of independent q is only 6 instead of 10,000. Then each 
q would be treated as if it were a separate coupling variable in Y*, Y A . 


4.2, Error Analysis of Fitted Response Surfaces 

For the current research, two tests were conducted to verify the accuracy of the 
fitted response surface. The first test involved the trivial case of fitting a quadratic 
response surface to a quadratic function. The second test involved the significantly more 
complex generalized Rosenbrock (Banana) function, 

NX 

J ( X ) = ^[lOO(x i+ , - x]y + (.Xj — !) ]• Equation 15 

;=i 

In both test cases, the response surfaces were generated using a random point placement 
scheme. 

A simple error analysis was performed to quantify differences in the fitted 
response surface model and the actual functions. For both functions, response surface 
function evaluations were compared to actual function evaluations using the following: 

Normalized Error = (F actua rF rs)/F actual- Equation 16 

The results of this error analysis are reported in Figure 13 through Figure 16, where NX 
is the number of design variables, and NS is the number of points required to fit the 
quadratic. As expected, the error associated with approximating a quadratic function 
with another quadratic function is smaller than the error found in the approximation to the 
Banana function. 
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Figure 13. Error analysis for NX=5, NS=21 
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Figure 14. Error analysis for NX=1 0, NS=66 
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Figure 15. Error analysis for NX=15, NS ::: 136 
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Figure 16. Error analysis for NX=20, NS=231 
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5. Design Point Placement 


A number of design point placement techniques were examined to determine the 
most efficient way to explore the design space. Since the BB analyses were somewhat 
computationally costly (especially in the case of structures and aerodynamics), and since 
there were a high number of input variables to some BB’s, traditional DOE patterns such 
as central composite or Box-Behnken were impractical. This led to the consideration and 
testing of some alternative point placement schemes: Hypersphere, Monte Carlo 
(random) and D-optimal Point Placement schemes. In accordance with DOE literature, 
the term design used in a DOE context refers to a particular point dispersion technique, 
and should be differentiated from design in the physical sense. 

5.1. Hypersphere Point Placement 

The first DOE scheme investigated involved points placed uniformly on a 
hypersphere in n-dimensions, coupled with an anti-bunching mechanism. The 
hypersphere placement of points reduces the likelihood of point bunching, a common 
phenomenon (especially in the comer locations) of a hypercube design. Additionally, 
since the nature of the experiments performed was deterministic, the number of points 
placed on the hypersphere surface was significantly lower than that of a CCD or Box- 
Behnken design. In fact, the number of points placed using the hypersphere is governed 
only by the minimum number of points necessary for a second order surface fit to the 
data. 

The method used for generating uniformly distributed random points on the 
surface of the hypersphere was adapted from [15]. To protect agai nst the possibility of a 
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point falling very closely to or on top of its nearest neighbor, an “anti-bunching” 
mechanism was developed. Overlapping points would in effect act like a single point, 
therefore reducing the total number of points by one. Since the goal here is to generate 
the fewest points required to fit a quadratic surface, an overlapping set of points would 
drop the number of points below the minimum required for a quadratic response surface 
fit. Additionally, the accuracy of the fitted surface is reduced when the points are 
bunched together. A complete description of the hypersphere point placement with 
antibunching algorithm can be found in Appendix A. 

Once implemented, the hypersphere point placement method appeared to exhibit 
unacceptable biasing toward the center of the design space, especially for models of 
higher dimension [16], This can be explained by examining Figure 17 through Figure 19. 



Dimension 


Figure 17. Volume of N-Dimeiisioual Hypercube [16] 
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Figure 18. Volume of N-Dimensional Hypersphere [16] 



Dimension 


Figure 19. Ratio of Hypershpere Volume to Hypercube Volume [16] 

Contrary to a cube in N-dimensions, a sphere does not continually increase in 
volume with increasing dimensionality. In fact, after as little as N = 5 the volume of the 
N-dimensional hypersphere actually decreases, resulting in the ratio of hypersphere 
volume to hypercube volume quickly approaches zero. One can safely deduce that the 
volume of the hypersphere at high dimensionality does not effectively overlap a 
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comparable hypercube volume. Thus, the design space of interest is poorly approximated 
by the hypersphere point placement scheme. The hypersphere point placement scheme 
was ultimately abandoned for these shortcomings. A detailed formulation for the volume 
of a hypersphere can be found in Appendix A. 

5.2. Random Point Placement 

The next point placement scheme involved a random Monte Carlo generation of 
input variables to each BB. The Matlab random number generator was used to generate a 
random vector §, whose elements lie between 0 and 1. Once mapped into [-1 1] 
coordinates, the random vectors of inputs were used as parameters for local optimization 
and formation of the response surface. 

The random point placement was extremely simple to implement and showed 
good dispersion of the input parameters. However, once the response surfaces were fit 
and examined, there was evidence that these quadratic models gave a poor representation 
of the actual black box, especially at the edges of the design space. The origin of this 
disparity could come from two possible sources. First, the minimum number of points 
required to fit a quadratic may have simply yielded a poor fit. Secondly, the extremities 
of the response surfaces appeared not to “anchor” themselves at the end bounds of the 
data. Of course, it is impossible to plot an 11 th dimensional surface, so these hypotheses 
are based on slicing the design space and viewing contour plots of various input 
combinations. 

For the reasons stated above, efforts were made to investigate the effects of using 
more than the required minimum number of points for a quadratic fit. Indeed, by 
increasing the number of points, the correlation of the RS surrogate to the actual BB 
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response was improved. However, it was found through experience that reducing the 
intervals over which RS’s were generated yielded better results. 

5.3. D-Optimal Point Placement 

To address the shortcomings of the random point placement scheme, a D-Optimal 
method was considered. The D-optimal design belongs to a class of computer generated 
designs first formulated in the 1970s. A D-optimal design is one that maximizes the 
determinant of X 1 X. It can be shown that det(X T X) is inversely proportional to the square 
of the volume of the confidence ellipsoid of the regression estimates of the linear model 
[12]. The volume of this confidence region is important because it reflects how well the 
set of coefficients are estimated. 

The cordexch (coordinate exchange) function in Matlab’s statistics toolbox was 
chosen to generate the D-optimal design. This is an iterative algorithm that operates by 
improving a starting design by making incremental changes to its elements [17]. The 
increments are the individual elements of the design matrix, cordexch requires the 
user to specify the number of inputs, the number of runs or design points, and the order of 
the model. An important feature of this type of design is that the user can specify the 
number of points placed, so long as that number is large enough to fit a second-order 
response surface through the data. 

Baseline response surfaces are critical in subsequent BLISS iterations because the 
design space is reduced, then focused on the region that contains the system optimum. 
Therefore, if a reduced interval does not contain the system optimum, poor results will 
encompass all other iterations. Ultimately, the baseline design that was chosen was a D- 
optimal design. This choice seemed most appropriate since it efficiently searched the 
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design space with a limited number of design points. The D-optimal design also seemed 
to provide reasonable consideration to the edges of the design space - an attribute 
common to CCD or Box-Behnken designs, but made possible with far fewer design 
points. 

Since BLISS-RS calls for the response surfaces to be fit through locally optimum 
data, it was imperative that all local optimizations were successfully converged. 
Therefore, a provision was implemented that allowed for the placement of a random 
design point whenever a particular D-Optimal design point failed to converge 
successfully. Such random points were generated until one was found that successfully 
converged. Thus, the D-Optimal design was ultimately augmented by the placement of 
randomly dispersed points 

5.4. Coded and Natural Variables 

A coordinate mapping is used to relate the coded variables generated by the point 
placement algorithm (ranging from -1 to 1) to the natural values of the variables in the 
physical problem. Figure 20 shows how such a mapping can be used to specify the 
natural coordinates in terms of their coded counterparts. 
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Figure 20. Linear mapping for Coded and Natural Variables 
In equation 17, §■ is the coded coordinate of a design point while X is the actual 
value. Notice that the high and low values for natural variables must be specified 
beforehand. These values may depend on the physical constraints of the problem or may 
simply be the region of interest for the experiment (e.g. the exploration interval). 


M-% l A - V 
~ §L - X L 


Equation 17 


6. Interval Reduction 

Once a system optimization is conducted, it is important to check the accuracy of 
the response surface fit. A simple calculation was used to compare values taken from the 
response surfaces to those produced directly from the BB analyses. First, the values {Z j 
Y j w} * are found by the system optimization. Then these values were used as BB input 
parameters, and the BB’s were re-optimized. The difference Ybb*~Yrs* was taken and 
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normalized by Ybb* to give a measure of the “goodness of fit” obtained in the quadratic 
response surface model. 

The measure of error due to lack of fit was used as a termination criteria. If, at the 
end of one iteration, the values Yrs* and Ybb* do not agree by a predetermined tolerance, 
the response surfaces are reformulated by reducing the intervals around system optimal 
values. 

Reducing the interval sizes has the general property of improving the fit of the 
response surface. Indeed, as the intervals are reduced, the fit becomes more and more 
linear. It is desired that the new interval be centered on the anticipated optimum value, 
and that the new interval makes physical sense. Since the Y and w variables are not 
subject to any clear physical limitations, their interval bounds are allowed to “drift” 
across the original interval bounds. However, the Z variables are subject to such physical 
limitations so their intervals must always lie within the original interval bounds. These 
differences lead to the use of two different i nterval reduction methods. 

Intervals on system level variables Z are reduced in such a way as to ensure the 
subsequent interval would fall within the physical bounds described by the original 
interval. For example, let XL; and XH; denote the respective lower and upper limits on 
the interval used to generate a response surface. If the termination tolerance is not 
satisfied, then the interval is reduced by 20% and centered on the system optimal value. 
If any of the new XL;+iis less than XL; , the value is returned to XL;. Similarly, if any of 
the XHj+i is larger than XH; the interval limit is returned to XH;. This operation shrinks 
the interval over which points are placed to generate the response surfaces. It should be 
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noted, however, that reducing the intervals in this manner does not ensure that the new 
interval is centered on the optimum found in the previous iteration. 

The factor by which the interval is reduced each iteration was chosen somewhat 
arbitrarily. Since the interval is reduced each iteration (without the possibility of being 
increased), it is important to reduce the interval size gradually to prevent premature 
convergence to an unreliable optimum 

The intervals for coupling variables Y and weight factors w were not subject to 
the limitation of falling within the original interval. Therefore, the intervals are allowed 
to shift solely on the basis of where the optimum was found on the previous iteration. 
This carries the added benefit of centering the RS’s in the middle of the design space with 
respect to Y and w, so satisfaction of the coupling constraints is always possible. 

A simple algorithm was adapted from [16] for re-sizing the intervals on Z and w. 
First, a factor K is defined for reducing the interval size. This factor is not directly 
applied to the current interval. Instead, an adjustment factor A , is determined by linear 
interpolation as: 


I C'-Opt, 

Equation 

C, ~ T 


where C, is the center of the current interval, Opt; is the optimum found in the current 
interval, and L; is the lower value of the current interval. The actual reduction factor can 
then be found as: 


R ( = A + (1 — A t )K Equation 19 

The result of this process is that the interval size will not be reduced if the optimum point 

is located on a boundary, while it will be reduced by a factor of K when the optimum 

point is located at the center. 
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Additionally, a history of optimum location with respect to bounds is kept. If the 
best design point hits the same bound twice, the interval is expanded by a fixed amount. 
If the optimum oscillates from upper to lower bounds between iterations, the interval is 
reduced by a fixed amount. 


7. BLISS-RS Results 

BLISS-RS was tested on the design of a supersonic business jet from a 1995 
AIAA Graduate Student Design Competition. The baseline model is shown in Figure 21 
with corresponding dimensions listed in Table 1. 



Figure 21. Baseline Model of’ Supersonic Business Jet 


Table 1. Baseline Geometrical Variables 




Initial Low Value 

Initial Value 

Initial High Value 


thickness ratio 

t/c 

.01 

0.075 

.1 

-- 

1 


3.5 

4 

6.5 

-- 



40 

45 

70 

degrees 

wing surface area 

Sw 

200 

400 

800 

ft 2 



40 

120 

125 

ft 2 

tail aspect ratio 

Mali 

3.5 

4.5 

5.5 

-- 

taper ratio 

A 

.1 

0.2 

.4 

— 
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For each BB, equation 14 was used to find the minimum number of points required for a 


quadratic RS. Table 2 lists the variables Z, Y*, and w that were input to each BB, and 
the corresponding number of design points placed. The weighting factors are labeled in 
such a way as to denote the BB number immediately following w, and the output number 
in parentheses. The weighting factor of one of the outputs of each BB was held constant 
(except for the performance BB since there is only one output). This imparts relative 
importance to the other weighting factors while also reducing the total number of 
variables. These were chosen as wl(l), the weighting factor for Wt, w2(2), the 
weighting factor for drag, and w3(2), the weighting factor for We, since it was obvious 
that minimizing these terms would have a beneficial effect on range. 


Table 2. RS Generation Data 


n 


Input V ariables 

Number 

of 

Inputs 

Number of 
Points (NS) 

Z 

Y* 

W 


1. W T 

2. W F 

3. 0 

t/c, ARw, Aw, 
Sri . Sht» 
ARht, 2. 

L, W E 

wl(2) 

wl(3) 


D 

2. Aero 

4. L 

5. D 

6. L/D 

t/c, h, M, 

A w , Sref, Sht, 
ARht, 7. 



14 

120 

3. Power 

1. SCF 

2. W E 

3. ESF 

h, M 

D 

w3(l) 

w3(3) 

5 

21 

4. Performance 

1 . Range 

h, M 

W T , 

W f , 

L/D, 

SFC 


6 

28 


For simplicity, the procedure was terminated after ten iterations. It was decided 
from previous executions that this number of iterations sufficiently demonstrated the 
convergence characteri stics of the system without excessive run time. 
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The total labor for each BB in Table 2 is NS times the number of iterations for 
convergence, the latter being 10 for the current research. If a concurrent processing 
environment were used, the elapsed time for each BB would be equal to the time for one 
NS generation. Therefore, the total elapsed time to generate all RS data would be equal 
to the elapsed time for the longest single BB optimization, and would scale linearly in 
proportion to the above longest time. 

Figure 22 through Figure 30 show the location of the optimum value for each of 
the system-level variables over the 1 0 iterations. Also shown are the intervals over which 
the response surfaces were formed per iteration. It is clear that by forcing all subsequent 
intervals to remain inside the original interval, some of the (i+l) th intervals are not 
centered on the i th optima. 
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Figure 22. Optimization History for t/c 
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Figure 23. Optimization History for Altitude 
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Figure 24. Optimization History for Mach Number 
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Figure 25. Optimization History for Wing Aspect Ratio 
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Figure 26. Optimization History for Wing Sweep (deg) 
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Figure 27. Optimization History for Wing Area (ft 2 ) 
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Figure 28. Optimization History for Tail Surface Area 
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Figure 29. Optimization History for Tail Aspect Ratio 
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Figure 30. Optimization History of Taper Ratio 
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Similarly, Figure 31 through Figure 39 show the location of the optimum value 
for each of the coupling variables, along with interval size for 10 iterations. In this case, 
the intervals are allowed to “drift” across the original interval bounds, and allows the 
(i+l) th interval to be centered on the i th optimum in all cases. This is a helpful 
characteristic when it comes to satisfying the coupling constraints (c}=0. The figures 
also demonstrate the workings of the algorithm given in equations 18 and 19. Intervals 
are reduced more aggressively if the optimum falls at the center of the interv al (a factor 
of 3 was used here). The intervals were expanded by a factor of two if the optimum 
landed at a boundary for two consecutive iterations, and the interval was reduced by a 
factor of two if oscillation occurred. 
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Figure 31. Optimization History for Total Aircraft Weight 
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Figure 32. Optimization History for Fuel Weight 
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Figure 33. Optimization History of “Twist” 
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Figure 34. Optimization History for Lift 
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Figure 35. Optimization History for Drag 
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Figure 37. Optimization History for SFC 
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Figure 38. Optimization History for Engine Weight 
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Figure 39. Optimization History for ESF 

The weighting factors are all plotted in Figure 40 for direct comparison. The 
values that the weighting factors converge to reflect the influence they have on the 
system objective. For instance, the weighting factor for lift, w2(l), is about -1, 
indicating that lift is being maximized (for the standard that positive numbers are 
minimized). Likewise, the weighting factor for SFC, w3(l), is 0.5, meaning SFC is being 
minimized. Both results are intuitively correct. 

Figure 40 also shows that the weighting factors for fuel weight, wl(2), and for 
ESF, w3(3), have not completely converged. This is clearly a product of the termination 


criteria chosen. 




~B~ wl (2) 
-A-~w1(3) 
-s«-~w2(1) 
-S~w2(3) 
■■■#■■■ w3(1) 
— i — w3(3) 


Figure 40. Optimization History for Weighting Factors, w 


The optimized ranges found by the RS representations and BB’s for each iteration 
are shown in Figure 41. In general, the response surface values closely resemble those 
obtained from the actual black box. The range of an actual aircraft would be substantially 
lower than those in shown in Figure 41. This is expected since constraints such as 
takeoff run length, climb rate requirements, engine-out conditions, and several aeroelastic 
limitations were not considered. 
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In order to assess the validity of the response surface models, the system level variables 
obtained in system optimization were used to re-optimize each BB. The corresponding 
outputs Ybb from this operation were compared to the Yrs obtained from the system 
optimizer. This comparison can be seen in Figure 42. 


50 






-B-Wtotal 
Wfuel 
-*-0 
-X-L 
-9K— D 
— © — L/D 
— +— SFC 
-m- Weng 
ESF 


Figure 42. Plot of Residuals for Coupling Constraints 

The data in Figure 42 were all normalized by the BB value in order to make a more direct 
comparison for the various coupling variables. As this graph show's, the interval 
shrinking technique does (in general) improve the accuracy of the response surface 
model. However, use of this method does not insure that each subsequent iteration will 
yield a better fit than the one before it. 

It is important to point out here that the comparison of the RS values to the BB 
values in Figure 42 is somewhat limited in scope. Ideally, the local variables should be 
represented by response surfaces of their own in order to show their behavior over the 
specified input intervals. One could then pull the information regarding local variables at 
the system optimum from these RS’s instead of re-optimizing local ly. Unfortunately, this 
would require a very large amount of data to be carried along through the B LISS process, 
and was therefore not performed. 
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Figure 43 through Figure 45 show the optimization histories for some local variables for 
structures, aerodynamics, and power black boxes, respectively. In Figure 43, the skin 
caliper thickness t s i through t S 3 are defined via Figure 6. The thickness are labeled 
“inner” to denote the spanwise location of the panels (inboard, or closest to the fuselage 
in this case). This figure also indicates, the values for t s i (inner) and t S 3 (inner) overlap. 
Additionally, the values for all three local variables shown do not deviate much from the 
initial values. Both anomalies can be attributed to the relatively relaxed tennination 
criteria chosen for the structures BB. This was done to speed up the construction of 
response surfaces. 



Figure 43. iteration Histories for Typical Structures BB Local Variables 
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Figure 45. Iteration History for Power BB Local Variable Thrust 
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8. Black Box Fidelity Improvements 


Since BLISS is modular by nature, the analyses included within each of the black 
boxes may be readily replaced. Indeed, the number of black boxes is also problem 
dependent. It is in this capacity that considerable efforts have been undertaken to 
demonstrate such capabilities. The following sections describe auxiliary work gone into 
BLISS-RS to improve the fidelity of submodule analyses. However, it should be noted 
that the results presented earlier do not pertain to these improvements. 

The most prominent improvement is the inclusion of a computational fluid 
dynamics package for aerodynamic analysis. This improved aerodynamics BB facilitates 
changes within the inner workings of the structures BB. The power and performance 
black boxes remained exactly the same as in the previous version of BLISS, aside from 
some trivial input/output modifications. 

Improved Aerodynamics Black Box 
The fidelity of the aerodynamics black box was improved by employing a 
computational fluid dynamics (CFD) code, coupled with grid-morphing software. Since 
the RS methodology calls for an optimization of local variables at each of the proposed 
design points and optimization of any local variables w'ould be extremely expensive 
computationally, it was decided that all local variables in the aerodynamics BB be 
eliminated. Thus the new aerodynamics BB is purely an analytical tool for finding the 
desired aerodynamic loads, namely Lift and Drag. 

The CFD code chosen for this task was CFL3D version 6.0, short for 
Computational Fluids Laboratory 3-Dimensional flow solver. CFL3D was originally 
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developed at by the Computational Fluids Laboratory at NASA Langley Research Center 


in the early 1980’s [ 19]. 

CFL3D solves the time-dependent conservation law form of the Reynolds- 
averaged Navier-Stokes equations. Some of the features of CFL3D are outlined below: 

® Semi-discrete finite volume approach for special discretization 
® Up-wind biasing for convective and pressure terms 

* Central differencing for shear stress and heat transfer terms 

* Implicit time advancement with ability to solve for steady or unsteady flows 

* Multigrid and mesh sequencing available for convergence acceleration 
® Multiple turbulence models available for 0, 1 ,or 2-equation models 

* Multiple-block toloplogies possible through 1-1 blocking, patching, overlapping, 
and embedding. 


CFL3D does not include any grid-generation software, so it was necessary to 
obtain the appropriate grid prior to execution. The baseline grid was generated at 
GEOLAB (Geometry Laboratory) at NASA Langely Research Center. The lab provides 
production and consultation services for computer aided design and numerical grid 
generation for various research fields. The baseline model chosen for BLISS-RS was that 
of the supersonic business jet used in previous BLISS work. The particular geometric 
parameters were obtained from the results of BLISS-98 optimization. 

Since each point on a RS represents a unique design, a unique surface and volume 
grid must be used for analysis at that point. Therefore, instead of reproducing the grid 
each time a new design was proposed, a grid-morphing program was used. The 
MASSOUD algorithm [20] chosen for this task. MASSOUD (Multidisciplinary 
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Aero/Struc Shape Optimization Using Deformation is a geometry parameterization tool 
that utilizes soft object animation algorithms used in computer graphics. It parameterizes 
the shape perturbations rather than the geometry itself, and relates grid deformation to 
aerodynamics shape design variables such as thickness, camber, twist, shear, and 
planform. The morphing capabilities available through MASSOUD are independent of 
grid topology, making it suitable for a variety of analysis codes such as CFD and CSM. 
Sensitivity derivatives are available in this software, and can be used for gradient-guided 
optimization. This algorithm is suitable for both low-fidelity (e.g., linear aerodynamics 
and equivalent laminated plate structures) and high-fidelity analysis tools (e.g., nonlinear 
CFD and detailed CSM modeling) [20]. A schematic of the data flow in this BB can be 
seen in Figure 46 below. 
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Figure 46. System Variable Flow for New BLISS-RS Aero BB 


Figure 46 also shows the output values for Lht and Lw, lift of the horizontal tail 
and wing, respectively. These values are used in a trim constraint during system 
optimization. Previous versions of BLISS performed this operation within the 
aerodynamics BB. 

Improved Structures Black Box 

In the upgraded version of BLISS, the actual lift distribution is passed along to the 
structures BB from the Aero BB via CFD analysis. The six coefficients (3 are used to 
describe the quadratic representation of the lift distribution. This is an improvement to 
the previous assumed elliptical /trapezoidal (Shrenk Approximation) distribution, thus the 
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accuracy of this distribution is greatly enhanced. Figure 47 shows the system variable 
flow for the improved structures BB. The notion of twist used in previous discussions 
has been abandoned here in favor of dA and dB, the vertical displacement at the wing tip 
leading edge and trailing edge, respectively. 



Figure 47. System Variable Flow for New BLISS-RS s ti nctures BB 


Another important modification to the structures BB is the implementation of an 
anti-flutter constraint during local optimization. The initial formulation of the flutter 
constraint (explained in detail in Appendix B) is based on the ratio of deformed wing 
bending moment, Mdr, to the rigid wing root bending moment, Mrr. However, a higher- 
fidelity flutter calculation is possible using FLAPS. 
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The changes to the aerodynamics and structures modules have an effect on the 
overall data flow of the system. Figure 48 shows the new data flow associated with 
implementation of the new aerodynamics and structures BB’s. Likewise, Figure 49 
diagrams the system optimization process associated with the new BB’s. 




Figure 48. New BLISS-RS System Data Flow 
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Figure 49. New BLISS RS System Optimizer 


9. Conclusions and Future Work 

The results presented herein reflect a change in the overall BLISS process: a shift 
from gradient-guided optimization methods to RSM based schemes. BL1SS-RS allows 
autonomous, distributed, and concurrent optimization in subsystems (disciplines) off-line. 
The method decouples local variables and constraints from the system level variables and 
coupling constraints, and optimizes the system for its objective and satisfies couplings. 
The response surface methodology was successfully implemented and tested on a 
conceptual aircraft design. By demonstrating the ability to use response surfaces in the 
BLISS-RS procedure, a number of possible improvements and follow-up work have 
emerged. 

First, the logical next step in the evolution of BLISS would be study its 
performance in a parallel-processing environment. The real gains to be made through 
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response surface methods are offset by the computational labor involved to create the 
RS’s. Since the RS’s generated in the current research were done serially, these gains 
were difficult to assess. 

The modular nature of BLISS-RS easily allows for improvements to the analyses 
local to each BB. Indeed, improvement in the fidelity of the black boxes was 
demonstrated in BLISS-98 [1-4], and again in BLISS-RS with the improved 
aerodynamics BB. A possible extension of this idea would be to increase the number of 
BB’s, therefore potentially increasing the number of all levels of design variables. 
Prospective new BB’s could include a stability and control analysis or even cost analysis. 
The process could also be adapted to reflect a more realistic design, for example by 
considering the entire flight envelope, or by more specific modeling. Therefore, BLISS- 
RS could be used in more advanced stages of the design process. 

Throughout the course of the research, many trial and error based decisions were 
conducted so there remain many possibilities for further investigation into “the road not 
taken”. The choice of point placement schemes was based primarily on judgmental 
criteria. Therefore, a definitive study should be conducted to find the best scheme for 
accurate results, and one that is practical to implement. The predominant case in this 
regard is the manner in which response surfaces are created from one iteration to the next. 
The BLISS-RS process can be made more efficient by using data points from previous 
iterations circumventing the need to produce a full set of points each iteration. This 
problem was addressed in the current research, but the dimensionality of the problem 
limited the number of points that could be retained, while still giving a set of points that 
fell within the interval of the next iteration. 
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Lastly, the user interface of BLISS is a prime candidate for impro vement. As the 
program stands, user defined inputs are inserted directly into the source code. This 
presents many challenges (experienced firsthand by the author) when trying to make 
small changes in baseline aircraft configuration, BB analyses, termination tolerances, etc. 
It will be essential that the interface becomes more user-friendly if BLISS-RS is to be 
used in actual design applications. 
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Appendix A. Hypersphere Point Placement 


This method generates uniformly distributed random points on the surface of the 
hypersphere, then employs an antibunching mechanism. A brief synopsis of the 
algorithm follows from [11] and [15]: 

1) . Generate normally distributed deviates X 1 ,X 2 ,...,X n with mean zero and variance 
one. 

2) . Obtain coordinates of points on the surface of the sphere by dividing each deviate 
by the root sum of the squares of the deviates: 


Xj X 2 


where r = ^jx; + X: + ...+ X;, 


Equation A1 


\r r r j 

Since the distribution function for the point (X,,X 2 ,...,X n ) has a density that only depends 


on the distance from the origin, it has uniform distribution when projected onto the 
surface of the sphere. 

To protect against the possibility of a point falling very closely to or on top of its 
nearest neighbor, an “anti-bunching” mechanism was developed. Overlapping points 
would in effect act like a single point, therefore reducing the total number of points by 
one. Since the goal here is to generate the fewest points required to fit a quadratic 
surface, an overlapping set of points would drop the number of points below the 
minimum (NS). Additionally, the accuracy of the fitted surface is reduced when the 
points are bunched together. 


A-l 



The solution for this problem is achieved by checking the distance of a point 
generated to its nearest neighbor. The distance between two points, A and B can be 
shown as 


dist AB 



(XA i - XB : ) 


Equation A2 


where XA and XB are the vectors containing the NX-coordi nates of points A and B, 
respectively. If the points falls too close to previous points, it is rejected and another 
point is generated until one meets a certain criteria. 

To introduce the notion of an anti-bunching criterion, suppose an N-dimensional 
coordinate system was partitioned into intervals (lower bound through upper bound) on 
each axis. After normalizing the lengths of the intervals to one, each is then divided into 
D divisions of equal length, 


S = — Equation A3 

D 

The quantity S can be thought of as a measure of distance from a point to its nearest 
neighbor expressed in units of the normalized interval. It may also be regarded as a 
measure of the mesh density. If each division on every axis were sampled at its center, 
the number of points generated would be NP = D v . After inverting this expressi on, the 
corresponding mesh density is 

S = NP N . Equation A4 

To prevent bunching, the distance between points can be limited to a l/li fraction of S, or 

S 

dlSt AB > — Equation A5 

h 
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Where h is a user-defined input, the anti-bunching factor. From this definition, a high 


anti-bunching factor will result in more bunching of points whereas lowering the value of 
h makes the point dispersion approach an even distribution, as shown in figures A1 
through A4. Also note, a single center point has been generated at the origin in each 


case. 



Figure AJ. Distribution of 50 points on a two-dimensional hypersphere with h=100. 
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Figure A2. Distribution of 50 points on a two-dimensional hypersphere with h=10. 



Figure A3. Distribution of 50 points on a two-dimensional hvpersphere with h=5. 



n 


Figure A4. Distribution of 50 points on a two-dimensional hypersphere with h=1.5. 

As expected, a lower anti-bunching factor results in a more uniform distribution. It is 
important to note that the anti-bunching factor must not fall below the value that gives a 
truly even distribution. If this were the case, all subsequent points generated would 
violate the distance criterion of equation A5. For the case set of 50 points given above 
(49 on the surface and one center run), this critical value is h cr it“l.l 12. 

The volume of a hypersphere can be found as from equation A6: 


Volume = 


it' 


r-+i 


where n is the number of independent variables, and 

Oj 

r(n) = f x n - l e~ x dx 

o 

or 


Equation A6 


Equation A7 
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r(u) 


T(n) 


n 



if n even 



if n odd. 


Equation A8 


Equation A9 
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Appendix B. Empirical Flutter Constraint Formulation 


The following formulation of an empirical flutter constraint was adapted from 
[11], The method limits excessive flutter by comparing the root bending moment of a 
twisted wing to that of a rigid wing. A more realistic formulation that includes vibration 
frequency of the wing is also included. 

To introduce the notion of a twist constraint, consider the change in local angle of 
attack, Aa, of a spanwise station along the wing: 

( dA - clB) 

A a = Equation B1 

C AB 

where dA and dB are the z-displacements of points A and B (leading and trailing edge 
locations, respectively) and c AB is the chord length measured from point A to point B. 

The aerodynamics module uses A a to modify the wing shape and to compute the 
resulting changes in aerodynamic loads. Suppose now that the pressure distribution 
along the span of the wing is given, either by assuming a distribution (as was done in 
BLISS 98) or by computing the distribution via CFD. The total pressure distribution is 
given by 

F av (x) = cr f (x) Equation B2 

where f(x) is the pressure distribution, and a is a constant that controls the lift ma g nitude. 
However, the lift given by equation B2 is valid only for rigid wings. Therefore, we must 
account for wing bending and twisting so that a chord at location x rotates per equation 
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B 3 . This rotation results in altering the local angle of attack. The effect of A a on total 
lift is 


AL(x) = a • fix) 


A a{x) 


Equation B3 


where Aa is given for the local chord in equation B3 and a,- is the angle of attack of the 
root chord. Now the total lift can be determined as an integrated sum of the pressure 
distribution and the change in lift: 


L = a ■ f f(x)C(x) + ^dx 

' «r r 


Equation B4 


where C(x) represents the wing geometry (root chord, tip chord, and span for a 
trapezoidal wing). In equation B4, the integration extends from root to tip. This 
integrated pressure must match the required lift Lo input to the structures BB from the 
system optimizer so that 



A trade off between structural weight and drag occurs because the wing loses 
some angle of attack to Aa. The loss in angle of attack intensifies outboard and the 
related loss in lift (AL in equation B3) results in a decrease in wing bending moment 
since the center of lift moves inboard. The reduction of wing bending moment is an 
advantage in terms of structural weight, but results in increased aerodynamic drag since 
the lift distribution departs from the minimum drag, elliptical distribution. However, for 
a transport aircraft, the drag penalty can be ignored if the wing is built to a jig shape that 
incorporates the negative twist anticipated at 1 g flight. In 1 g flight, the wing flexes to 
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an ideal shape as if it were rigid, but still relieves the bending load in a pull-up maneuver. 
Hence, we can exploit the wing flexing as a natural mechanism to reduce the pull-up 
bending and thus save weight to improve range. 

The deformed wing root bending moment can be defined as 

M dr = a ' j~ fi x ) ’ 

Similarly, a rigid wing root bending moment can be written as 

M rr = a ■ J f(x)- C(x)xdx Equation B7 

Since Aa <0 for an aft-swept wing, r can be introduced as 

r = —cl Equation B8 

m rr 

When the value of r is less than unity, the wing is relieved of the bending load. For local 
optimization, r is used as a constraint to prevent excessive flutter. For example, the 
simplest form of the constraint would be a bound on r: 

r < r Equation B9 


C(x) + 


A a 


a „ 


xdx 


Equation B6 


where r~ could be a values such as 0.9, it which case 10% bending relief would be 
allowed. A better approximation can be formulated based on the formulas of vibration 
frequency. The formulas for bending and torsional frequency are given in equations BIO 
and B 1 1 , respectively. 


(O 


b 



Equation BIO 
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k , 

2 \7 


Equation Bll 


where Pi and P 2 are constants of proportionality, kb is the bending stiffness corresponding 
to the wing tip unit torque, k t is the torsion stiffness corresponding to the wing tip unit 
torque, m is wing structural mass, and J is the wing moment of inertia associated with 
rotation about the elastic axis. J can be expressed in terms of m and radius of gyration, R, 

by: 

J = mR Equation B12 

If the wing is free of flutter at some reference state (designated with a ~), limits can be set 
around (Ob and (Ot (ub and u t ) so that 

g = — - — (l— U,)<0 Equation B13 


g = (l— u,)-—<0 Equation B1 4 

These constraints keep 0 % within ± Ub of co 6 . The actual value of Ub is somewhat 
arbitrary, for example, 0.1. Similar constraints can be generated for (Ot. The reference 
state ~ can be detennined by executing BLISS with the constraint in equation B9. Then 
the constraint given in equation B9 may be replaced with equations B13 and B14. Note 
that in equations B13 and B14, m s , kb, and k t are the only items that need to be computed 
since Pi, P 2 , and R cancel out (under a simplifying assumption that they remain 
approximately constant). 
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Figure Cl. BLISS-RS Results for Local Variables 
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Figure C2. BLISS-RS Results for System Variables 
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Figure C3. BLISS-RS Results for Coupling Variables 
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engineering project involving a number of specialty groups that might be geographically dispersed, and it exploits the contemporary computing 
technology of massively concurrent and distributed processing. The report includes a demonstration test case of supersonic business jet design. 
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